package 'xlsx' successfully unpacked and MD5 sums checked
The downloaded binary packages are in
C:\Users\cvasquezv\AppData\Local\Temp\Rtmpy0PSlr\downloaded_packages
Diseño completamente al azar en arreglo de parcelas divididas
modelo.dca1 <-lm(rdto ~ A * B + rep/A - rep, data = data)modelo.dca2 <-lm(rdto ~ A * B + rep:A - rep, data = data)modelo.dca3 <-lm(rdto ~ A + B + rep/A - rep, data = data)modelo.dca4 <-lm(rdto ~ A + B + rep:A - rep, data = data)
Nota
La expresión “rep/A” es para considerar la interacción como efecto fijo.
La expresión “rep:A” es para considerar la interacción como efecto aleatorio.
broom::glance(modelo.dca1) %>%bind_rows(broom::glance(modelo.dca2), broom::glance(modelo.dca3), broom::glance(modelo.dca4)) %>% dplyr::mutate(Modelo =c("A * B + rep/A","A * B + rep:A","A + B + rep/A","A + B + rep:A")) %>% dplyr::select(Modelo, AIC, BIC) %>% dplyr::arrange(BIC) %>% dplyr::mutate(Mérito =1:n()) %>% dplyr::relocate(Mérito, Modelo) %>%gt()
Mérito
Modelo
AIC
BIC
1
A * B + rep/A
84.76335
94.55744
2
A * B + rep:A
84.76335
94.55744
3
A + B + rep/A
121.70347
129.71682
4
A + B + rep:A
121.70347
129.71682
Definición del modelo
modelo.dca <-lm(rdto ~ A * B + rep/A - rep, data = data)
Durbin-Watson test
data: modelo.dca
DW = 2.4588, p-value = 0.8145
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dca))
Shapiro-Wilk normality test
data: rstudent(modelo.dca)
W = 0.95347, p-value = 0.482
ad.test(rstudent(modelo.dca))
Anderson-Darling normality test
data: rstudent(modelo.dca)
A = 0.28753, p-value = 0.5783
lillie.test(rstudent(modelo.dca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dca)
D = 0.12222, p-value = 0.6786
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dca)
D = 0.12079, p-value = 0.9274
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dca))
Cramer-von Mises normality test
data: rstudent(modelo.dca)
W = 0.035902, p-value = 0.7384
pearson.test(rstudent(modelo.dca))
Pearson chi-square normality test
data: rstudent(modelo.dca)
P = 3, p-value = 0.5578
sf.test(rstudent(modelo.dca))
Shapiro-Francia normality test
data: rstudent(modelo.dca)
W = 0.94303, p-value = 0.2775
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 4.009865, Df = 1, p = 0.045235
bptest(modelo.dca)
studentized Breusch-Pagan test
data: modelo.dca
BP = 14.223, df = 9, p-value = 0.1146
bptest(modelo.dca, studentize = F)
Breusch-Pagan test
data: modelo.dca
BP = 10.037, df = 9, p-value = 0.3475
olsrr::ols_test_breusch_pagan(modelo.dca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
-----------------------------
DF = 1
Chi2 = 4.009865
Prob > Chi2 = 0.04523478
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dca %>%gvlma()
Call:
lm(formula = rdto ~ A * B + rep/A - rep, data = data)
Coefficients:
(Intercept) Aa2 Bb2 Bb3 Aa2:Bb2 Aa2:Bb3
2.6667 1.7778 -1.6667 7.6667 9.3333 -10.6667
Aa1:rep2 Aa2:rep2 Aa1:rep3 Aa2:rep3
1.3333 0.3333 1.6667 -0.6667
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 3.20107 0.5248 Assumptions acceptable.
Skewness 0.07893 0.7788 Assumptions acceptable.
Kurtosis 0.25987 0.6102 Assumptions acceptable.
Link Function 2.78083 0.0954 Assumptions acceptable.
Heteroscedasticity 0.08145 0.7753 Assumptions acceptable.
\(H_1: \text{En al menos una interacción entre el factor A y el factor B el } (\tau\beta) \text{ es diferente a los demás.}\)
\(H_1: (\tau\beta)_{ij} \neq 0\text{; en al menos una interacción entre el factor A y el factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
anova(modelo.dca, test ="F")
Analysis of Variance Table
Response: rdto
Df Sum Sq Mean Sq F value Pr(>F)
A 1 0.222 0.222 0.0516 0.8259782
B 2 29.778 14.889 3.4581 0.0827438 .
A:B 2 300.444 150.222 34.8903 0.0001119 ***
A:rep 4 6.222 1.556 0.3613 0.8296470
Residuals 8 34.444 4.306
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Importante
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(A:parcelagrande)” en el caso de efectos aleatorios o “Error(A:rep)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(rdto ~ A * B +Error(rep/A), data = data) -> aov.dcasummary(aov.dca)
Error: rep
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 2 2.111 1.056
Error: rep:A
Df Sum Sq Mean Sq F value Pr(>F)
A 1 0.222 0.2222 0.108 0.774
Residuals 2 4.111 2.0556
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
B 2 29.78 14.89 3.458 0.082744 .
A:B 2 300.44 150.22 34.890 0.000112 ***
Residuals 8 34.44 4.31
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
broom::tidy(aov.dca)
# A tibble: 6 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 rep Residuals 2 2.11 1.06 NA NA
2 rep:A A 1 0.222 0.222 0.108 0.774
3 rep:A Residuals 2 4.11 2.06 NA NA
4 Within B 2 29.8 14.9 3.46 0.0827
5 Within A:B 2 300. 150. 34.9 0.000112
6 Within Residuals 8 34.4 4.31 NA NA
Valor de la tabla de F para el factor A con una significancia de 0.05.
qf(0.95, 1, 2)
[1] 18.51282
Valor de la tabla de F para el factor B con una significancia de 0.05.
qf(0.95, 2, 8)
[1] 4.45897
Valor de la tabla de F para la interacción A:B con una significancia de 0.05.
Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor A tienen un efecto sobre el rendimiento estadísticamente similar a 0.
Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor B tienen un efecto sobre el rendimiento estadísticamente similar a 0.
Con respecto a la interacción entre el Factor A y Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos una interacción entre un nivel del factor A y un nivel del factor B existe un efecto de antagonismo o sinergismo sobre el rendimiento.
agricolae::cv.model(modelo.dca)
[1] 35.91317
Comparaciones de medias para los efectos principales del Factor A
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta A, # Cambiar según nombre de variable independienteDFerror =get_df_ea_spplot(aov.dca), MSerror =get_mse_ea_spplot(aov.dca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ A
LSD t Test for rdto
Mean Square Error: 2.055556
A, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
a1 5.666667 4.690416 9 3.610399 7.722934 1 15
a2 5.888889 4.935698 9 3.832621 7.945157 1 14
Alpha: 0.05 ; DF Error: 2
Critical Value of t: 4.302653
least Significant Difference: 2.908002
Treatments with the same letter are not significantly different.
rdto groups
a2 5.888889 a
a1 5.666667 a
Comparaciones de medias para las interacciones
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
B1 vs B3:
\(H_0: \mu_{B1} - \mu_{B3} = 0\)
\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)
B2 vs B3:
\(H_0: \mu_{B2} - \mu_{B3} = 0\)
\(H_1: \mu_{B2} - \mu_{B3} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
Análisis de varianza para interacción de dos factores con el paquete phia
Comparación de los niveles de B dentro de cada nivel de A
Comparaciones de medias de los niveles de B dentro de cada nivel de A
FA_A1 <- data %>%filter(A=="a1")FA_A2 <- data %>%filter(A=="a2")fcomp1 <-function(x){ comp <-LSD.test(x$rdto, # Cambiar según nombre de variable respuesta x$B, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dca), MSerror = dvmisc::get_mse(modelo.dca),alpha =0.05,group=TRUE,main =NULL,console=FALSE)return(comp[[5]] %>%rename("rdto"="x$rdto") # Cambiar según nombre de variable respuesta )}
Comparaciones de los niveles de B dentro del nivel A1
fcomp1(FA_A1)
rdto groups
b3 11.333333 a
b1 3.666667 b
b2 2.000000 b
Comparaciones de los niveles de B dentro del nivel A2
fcomp1(FA_A2)
rdto groups
b2 12.000000 a
b1 4.333333 b
b3 1.333333 b
Comparaciones de medias de los niveles de B dentro de cada nivel de A
FB_B1 <- data %>%filter(B=="b1")FB_B2 <- data %>%filter(B=="b2")FB_B3 <- data %>%filter(B=="b3")fcomp2 <-function(x){ comp <-LSD.test(x$rdto, # Cambiar según nombre de variable respuesta x$A, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dca), MSerror = dvmisc::get_mse(modelo.dca),alpha =0.05,group=TRUE,main =NULL,console=FALSE)return(comp[[5]] %>%rename("rdto"="x$rdto") # Cambiar según nombre de variable respuesta )}
Comparaciones de los niveles de A dentro del nivel B1
fcomp2(FB_B1)
rdto groups
a2 4.333333 a
a1 3.666667 a
Comparaciones de los niveles de A dentro del nivel B2
fcomp2(FB_B2)
rdto groups
a2 12 a
a1 2 b
Comparaciones de los niveles de A dentro del nivel B2
fcomp2(FB_B3)
rdto groups
a1 11.333333 a
a2 1.333333 b
DCA en arreglo factorial de parcelas divididas con el paquete ExpDes
------------------------------------------------------------------------
Legend:
FACTOR 1 (plot): Variedad
FACTOR 2 (split-plot): Densidad
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr>Fc
Variedad 1 0.22 0.222 0.143 0.72466
Error a 4 6.22 1.556
Densidad 2 29.78 14.889 3.458 0.08274
Variedad*Densidad 2 300.44 150.222 34.890 0.00011
Error b 8 34.44 4.306
Total 17 371.11
------------------------------------------------------------------------
CV 1 = 21.58649 %
CV 2 = 35.91317 %
------------------------------------------------------------------------
Shapiro-Wilk normality test (Error b)
p-value: 0.8685814
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
Significant interaction: analyzing the interaction
------------------------------------------------------------------------
Analyzing Variedad inside of each level of Densidad
------------------------------------------------------------------------
DF SS MS Fc p.value
Variedad : Densidad b1 1.00000 0.666667 0.666667 0.196721 0.666407
Variedad : Densidad b2 1.00000 150.000000 150.000000 44.262289 0.000046
Variedad : Densidad b3 1.00000 150.000000 150.000000 44.262289 0.000046
Pooled Error 10.46818 35.475483 3.388889 NA NA
------------------------------------------------------------------------
Variedad inside of Densidad b1
------------------------------------------------------------------------
According to F test, the means of this factor are not distinct.
------------------------------------------------------------------------
Levels Means
1 a1 3.666667
2 a2 4.333333
------------------------------------------------------------------------
Variedad inside of Densidad b2
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a a2 12
b a1 2
------------------------------------------------------------------------
Variedad inside of Densidad b3
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a a1 11.33333
b a2 1.333333
------------------------------------------------------------------------
Analyzing Densidad inside of each level of Variedad
------------------------------------------------------------------------
DF SS MS Fc p.value
Densidad : Variedad a1 2 148.66667 74.333333 17.26451 0.001252
Densidad : Variedad a2 2 181.55556 90.777778 21.08387 0.000647
Error b 8 34.44444 4.305555 NA NA
------------------------------------------------------------------------
Densidad inside of Variedad a1
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a b3 11.33333
b b1 3.666667
b b2 2
------------------------------------------------------------------------
------------------------------------------------------------------------
Densidad inside of Variedad a2
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a b2 12
b b1 4.333333
b b3 1.333333
------------------------------------------------------------------------
------------------------------------------------------------------------
Importante
El paquete ExpDes usa la expresión “rep:A”, es decir, considera la interacción como efecto aleatorio.
Diseño de bloques completos al azar en arreglo de parcelas divididas
modelo.dbca1 <-lm(rdto ~ A * B + bloque/A + bloque, data = data)modelo.dbca2 <-lm(rdto ~ A * B + bloque:A + bloque, data = data)modelo.dbca3 <-lm(rdto ~ A + B + bloque/A + bloque, data = data)modelo.dbca4 <-lm(rdto ~ A + B + bloque:A + bloque, data = data)
Durbin-Watson test
data: modelo.dbca
DW = 2.4588, p-value = 0.8145
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.95347, p-value = 0.482
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.28753, p-value = 0.5783
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.12222, p-value = 0.6786
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.12079, p-value = 0.9274
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.035902, p-value = 0.7384
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 3, p-value = 0.5578
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.94303, p-value = 0.2775
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 4.009865, Df = 1, p = 0.045235
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 14.223, df = 9, p-value = 0.1146
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 10.037, df = 9, p-value = 0.3475
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
-----------------------------
DF = 1
Chi2 = 4.009865
Prob > Chi2 = 0.04523478
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dbca %>%gvlma()
Call:
lm(formula = rdto ~ A * B + bloque/A + bloque, data = data)
Coefficients:
(Intercept) Aa2 Bb2 Bb3 bloque2 bloque3
2.667 1.778 -1.667 7.667 1.333 1.667
Aa2:Bb2 Aa2:Bb3 Aa2:bloque2 Aa2:bloque3
9.333 -10.667 -1.000 -2.333
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 3.20107 0.5248 Assumptions acceptable.
Skewness 0.07893 0.7788 Assumptions acceptable.
Kurtosis 0.25987 0.6102 Assumptions acceptable.
Link Function 2.78083 0.0954 Assumptions acceptable.
Heteroscedasticity 0.08145 0.7753 Assumptions acceptable.
\(H_1: \text{En al menos una interacción entre el factor A y el factor B el } (\tau\beta) \text{ es diferente a los demás.}\)
\(H_1: (\tau\beta)_{ij} \neq 0\text{; en al menos una interacción entre el factor A y el factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
anova(modelo.dbca, test ="F")
Analysis of Variance Table
Response: rdto
Df Sum Sq Mean Sq F value Pr(>F)
A 1 0.222 0.222 0.0516 0.8259782
B 2 29.778 14.889 3.4581 0.0827438 .
bloque 2 2.111 1.056 0.2452 0.7882486
A:B 2 300.444 150.222 34.8903 0.0001119 ***
A:bloque 2 4.111 2.056 0.4774 0.6369845
Residuals 8 34.444 4.306
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Importante
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(A:parcelagrande)” en el caso de efectos aleatorios o “Error(A:bloque)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(rdto ~ A * B + bloque +Error(bloque/A), data = data) -> aov.dbcasummary(aov.dbca)
Error: bloque
Df Sum Sq Mean Sq
bloque 2 2.111 1.056
Error: bloque:A
Df Sum Sq Mean Sq F value Pr(>F)
A 1 0.222 0.2222 0.108 0.774
Residuals 2 4.111 2.0556
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
B 2 29.78 14.89 3.458 0.082744 .
A:B 2 300.44 150.22 34.890 0.000112 ***
Residuals 8 34.44 4.31
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
data %>%with(agricolae::sp.plot(block = bloque,pplot = A,splot = B,Y = rdto))
ANALYSIS SPLIT PLOT: rdto
Class level information
A : a1 a2
B : b1 b2 b3
bloque : 1 2 3
Number of observations: 18
Analysis of Variance Table
Response: rdto
Df Sum Sq Mean Sq F value Pr(>F)
bloque 2 2.111 1.056 NaN NaN
A 1 0.222 0.222 0.1081 0.7735446
Ea 2 4.111 2.056 NaN NaN
B 2 29.778 14.889 3.4581 0.0827438 .
A:B 2 300.444 150.222 34.8903 0.0001119 ***
Eb 8 34.444 4.306 NaN NaN
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cv(a) = 24.8 %, cv(b) = 35.9 %, Mean = 5.777778
broom::tidy(aov.dbca)
# A tibble: 6 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 bloque bloque 2 2.11 1.06 NA NA
2 bloque:A A 1 0.222 0.222 0.108 0.774
3 bloque:A Residuals 2 4.11 2.06 NA NA
4 Within B 2 29.8 14.9 3.46 0.0827
5 Within A:B 2 300. 150. 34.9 0.000112
6 Within Residuals 8 34.4 4.31 NA NA
Valor de la tabla de F para el factor A con una significancia de 0.05.
qf(0.95, 1, 2)
[1] 18.51282
Valor de la tabla de F para el factor B con una significancia de 0.05.
qf(0.95, 2, 8)
[1] 4.45897
Valor de la tabla de F para la interacción A:B con una significancia de 0.05.
Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor A tienen un efecto sobre el rendimiento estadísticamente similar a 0.
Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor B tienen un efecto sobre el rendimiento estadísticamente similar a 0.
Con respecto a la interacción entre el Factor A y Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos una interacción entre un nivel del factor A y un nivel del factor B existe un efecto de antagonismo o sinergismo sobre el rendimiento.
agricolae::cv.model(modelo.dbca)
[1] 35.91317
Comparaciones de medias para los efectos principales del Factor A
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta A, # Cambiar según nombre de variable independienteDFerror =get_df_ea_spplot(aov.dbca), MSerror =get_mse_ea_spplot(aov.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ A
LSD t Test for rdto
Mean Square Error: 2.055556
A, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
a1 5.666667 4.690416 9 3.610399 7.722934 1 15
a2 5.888889 4.935698 9 3.832621 7.945157 1 14
Alpha: 0.05 ; DF Error: 2
Critical Value of t: 4.302653
least Significant Difference: 2.908002
Treatments with the same letter are not significantly different.
rdto groups
a2 5.888889 a
a1 5.666667 a
Comparaciones de medias para los efectos principales del Factor B
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta B, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ B
LSD t Test for rdto
Mean Square Error: 4.305556
B, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
b1 4.000000 1.414214 6 2.046565 5.953435 2 6
b2 7.000000 5.656854 6 5.046565 8.953435 1 14
b3 6.333333 5.853774 6 4.379898 8.286769 1 15
Alpha: 0.05 ; DF Error: 8
Critical Value of t: 2.306004
least Significant Difference: 2.762575
Treatments with the same letter are not significantly different.
rdto groups
b2 7.000000 a
b3 6.333333 ab
b1 4.000000 b
Precaución
Si el análisis de varianza arroja que los niveles de un factor son estadísticamente similares entre sí, entonces no es necesario realizar una prueba de comparación de medias y por ende todos estos niveles pertenecen a un mismo grupo de significancia “a”.
Comparaciones de medias para las interacciones
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
B1 vs B3:
\(H_0: \mu_{B1} - \mu_{B3} = 0\)
\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)
B2 vs B3:
\(H_0: \mu_{B2} - \mu_{B3} = 0\)
\(H_1: \mu_{B2} - \mu_{B3} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
Análisis de varianza para interacción de dos factores con el paquete phia
Comparación de los niveles de B dentro de cada nivel de A
Comparaciones de medias de los niveles de B dentro de cada nivel de A
FA_A1 <- data %>%filter(A=="a1")FA_A2 <- data %>%filter(A=="a2")fcomp1 <-function(x){ comp <-LSD.test(x$rdto, # Cambiar según nombre de variable respuesta x$B, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=FALSE)return(comp[[5]] %>%rename("rdto"="x$rdto") # Cambiar según nombre de variable respuesta )}
Comparaciones de los niveles de B dentro del nivel A1
fcomp1(FA_A1)
rdto groups
b3 11.333333 a
b1 3.666667 b
b2 2.000000 b
Comparaciones de los niveles de B dentro del nivel A2
fcomp1(FA_A2)
rdto groups
b2 12.000000 a
b1 4.333333 b
b3 1.333333 b
Comparaciones de medias de los niveles de B dentro de cada nivel de A
FB_B1 <- data %>%filter(B=="b1")FB_B2 <- data %>%filter(B=="b2")FB_B3 <- data %>%filter(B=="b3")fcomp2 <-function(x){ comp <-LSD.test(x$rdto, # Cambiar según nombre de variable respuesta x$A, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group=TRUE,main =NULL,console=FALSE)return(comp[[5]] %>%rename("rdto"="x$rdto") # Cambiar según nombre de variable respuesta )}
Comparaciones de los niveles de A dentro del nivel B1
fcomp2(FB_B1)
rdto groups
a2 4.333333 a
a1 3.666667 a
Comparaciones de los niveles de A dentro del nivel B2
fcomp2(FB_B2)
rdto groups
a2 12 a
a1 2 b
Comparaciones de los niveles de A dentro del nivel B3
fcomp2(FB_B3)
rdto groups
a1 11.333333 a
a2 1.333333 b
DBCA en arreglo factorial de parcelas divididas con el paquete ExpDes
------------------------------------------------------------------------
Legend:
FACTOR 1 (plot): Variedad
FACTOR 2 (split-plot): Densidad
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of Variance Table
------------------------------------------------------------------------
DF SS MS Fc Pr(>Fc)
Variedad 1 0.22 0.222 0.108 0.773545
Block 2 2.11 1.056 0.514 0.660714
Error a 2 4.11 2.056
Densidad 2 29.78 14.889 3.458 0.082744 .
Variedad*Densidad 2 300.44 150.222 34.890 0.000112 ***
Error b 8 34.44 4.306
Total 17 371.11
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------------------------
CV 1 = 24.8144 %
CV 2 = 35.91317 %
No significant interaction: analyzing the simple effects
------------------------------------------------------------------------
Variedad
According to F test, the means of this factor are not different.
------------------------------------------------------------------------
Levels Means
1 a1 5.666667
2 a2 5.888889
------------------------------------------------------------------------
Densidad
According to F test, the means of this factor are not different.
------------------------------------------------------------------------
Levels Means
1 b1 4.000000
2 b2 7.000000
3 b3 6.333333
------------------------------------------------------------------------
Significant interaction: analyzing the interaction
------------------------------------------------------------------------
Analyzing Variedad inside of each level of Densidad
------------------------------------------------------------------------
DF SS MS Fc p.value
Variedad : Densidad b1 1.000000 0.666667 0.666667 0.18750 0.674199
Variedad : Densidad b2 1.000000 150.000000 150.000000 42.18749 0.000069
Variedad : Densidad b3 1.000000 150.000000 150.000000 42.18749 0.000069
Pooled Error 9.996678 35.543748 3.555556 NA NA
------------------------------------------------------------------------
Variedad inside of Densidad b1
------------------------------------------------------------------------
According to F test, the means of this factor are not different.
------------------------------------------------------------------------
Levels Means
1 a1 3.666667
2 a2 4.333333
------------------------------------------------------------------------
Variedad inside of Densidad b2
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a a2 12
b a1 2
------------------------------------------------------------------------
Variedad inside of Densidad b3
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a a1 11.33333
b a2 1.333333
------------------------------------------------------------------------
Analyzing Densidad inside of each level of Variedad
------------------------------------------------------------------------
DF SS MS Fc p.value
Densidad : Variedad a1 2 148.66667 74.333333 17.26451 0.001252
Densidad : Variedad a2 2 181.55556 90.777778 21.08387 0.000647
Error b 8 34.44444 4.305555 NA NA
------------------------------------------------------------------------
Densidad inside of Variedad a1
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a b3 11.33333
b b1 3.666667
b b2 2
------------------------------------------------------------------------
------------------------------------------------------------------------
Densidad inside of Variedad a2
------------------------------------------------------------------------
Tukey's test
------------------------------------------------------------------------
Groups Treatments Means
a b2 12
b b1 4.333333
b b3 1.333333
------------------------------------------------------------------------
------------------------------------------------------------------------
Importante
El paquete ExpDes usa la expresión “bloque:A”, es decir, considera la interacción como efecto aleatorio.
Diseño cuadrado latino en arreglo de parcelas divididas
modelo.dcl1 <-lm(PESO ~ row + col + FA * FB + row/FA, data = data)modelo.dcl2 <-lm(PESO ~ row + col + FA + FB + row/FA, data = data)modelo.dcl3 <-lm(PESO ~ row + FA * FB + row/FA, data = data)modelo.dcl4 <-lm(PESO ~ row + FA + FB + row/FA, data = data)modelo.dcl5 <-lm(PESO ~ col + FA * FB + row/FA, data = data)modelo.dcl6 <-lm(PESO ~ col + FA + FB + row/FA, data = data)
Durbin-Watson test
data: modelo.dcl
DW = 2.75, p-value = 0.3847
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del peso son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del peso es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del peso es similar a la función normal}\)
shapiro.test(rstudent(modelo.dcl))
Shapiro-Wilk normality test
data: rstudent(modelo.dcl)
W = 0.97736, p-value = 0.9709
ad.test(rstudent(modelo.dcl))
Anderson-Darling normality test
data: rstudent(modelo.dcl)
A = 0.16739, p-value = 0.9154
lillie.test(rstudent(modelo.dcl))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dcl)
D = 0.1384, p-value = 0.7587
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dcl)
D = 0.17491, p-value = 0.7975
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dcl))
Cramer-von Mises normality test
data: rstudent(modelo.dcl)
W = 0.028196, p-value = 0.8558
pearson.test(rstudent(modelo.dcl))
Pearson chi-square normality test
data: rstudent(modelo.dcl)
P = 2, p-value = 0.5724
sf.test(rstudent(modelo.dcl))
Shapiro-Francia normality test
data: rstudent(modelo.dcl)
W = 0.98564, p-value = 0.9892
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del peso es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del peso es constante con respecto a los valores ajustados del peso
\(H_1\): La varianza del peso no es constante con respecto a los valores ajustados del peso
ncvTest(modelo.dcl)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.427099, Df = 1, p = 0.51342
bptest(modelo.dcl)
studentized Breusch-Pagan test
data: modelo.dcl
BP = 12, df = 7, p-value = 0.1006
bptest(modelo.dcl, studentize = F)
Breusch-Pagan test
data: modelo.dcl
BP = 3.1406, df = 7, p-value = 0.8717
olsrr::ols_test_breusch_pagan(modelo.dcl)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : PESO
Variables: fitted values of PESO
Test Summary
----------------------------
DF = 1
Chi2 = 0.427099
Prob > Chi2 = 0.5134159
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del peso es constante con respecto a los valores ajustados del peso.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al peso, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dcl %>%gvlma()
Call:
lm(formula = PESO ~ row + FA * FB + row/FA, data = data)
Coefficients:
(Intercept) row2 FAA2 FBB2 FBB3 FAA2:FBB2
1.017e+01 6.667e-01 6.833e+00 6.500e+00 9.000e+00 6.527e-15
FAA2:FBB3 row2:FAA2
-1.350e+01 -6.667e-01
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 4.702e+00 0.31921 Assumptions acceptable.
Skewness 2.395e-32 1.00000 Assumptions acceptable.
Kurtosis 1.090e+00 0.29644 Assumptions acceptable.
Link Function 3.547e+00 0.05964 Assumptions acceptable.
Heteroscedasticity 6.503e-02 0.79872 Assumptions acceptable.
\(H_1: \text{En al menos una interacción entre el factor A y el factor B el } (\tau\beta) \text{ es diferente a los demás.}\)
\(H_1: (\tau\beta)_{ij} \neq 0\text{; en al menos una interacción entre el factor A y el factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(A:parcelagrande)” en el caso de efectos aleatorios o “Error(A:bloque)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
aov(PESO ~ FA * FB + row +Error(row/FA), data = data) -> aov.dclsummary(aov.dcl)
Error: row
Df Sum Sq Mean Sq
row 1 0.3333 0.3333
Error: row:FA
Df Sum Sq Mean Sq F value Pr(>F)
FA 1 12.000 12.000 36 0.105
Residuals 1 0.333 0.333
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
FB 2 87.17 43.58 32.69 0.00332 **
FA:FB 2 121.50 60.75 45.56 0.00177 **
Residuals 4 5.33 1.33
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
broom::tidy(aov.dcl)
# A tibble: 6 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 row row 1 0.333 0.333 NA NA
2 row:FA FA 1 12.0 12.0 36.0 0.105
3 row:FA Residuals 1 0.333 0.333 NA NA
4 Within FB 2 87.2 43.6 32.7 0.00332
5 Within FA:FB 2 122. 60.8 45.6 0.00177
6 Within Residuals 4 5.33 1.33 NA NA
Valor de la tabla de F para el factor A con una significancia de 0.05.
qf(0.95, 1, 1)
[1] 161.4476
Valor de la tabla de F para el factor B con una significancia de 0.05.
qf(0.95, 2, 4)
[1] 6.944272
Valor de la tabla de F para la interacción A:B con una significancia de 0.05.
Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor A tienen un efecto sobre el rendimiento estadísticamente similar a 0.
Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tienen un efecto sobre el rendimiento estadísticamente diferente a 0.
Con respecto a la interacción entre el Factor A y Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos una interacción entre un nivel del factor A y un nivel del factor B existe un efecto de antagonismo o sinergismo sobre el rendimiento.
Comparaciones de medias para los efectos principales del Factor A
data %>%with(LSD.test( PESO, # Cambiar según nombre de variable respuesta FA, # Cambiar según nombre de variable independienteDFerror =get_df_ea_spplot(aov.dcl), MSerror =get_mse_ea_spplot(aov.dcl),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: PESO ~ FA
LSD t Test for PESO
Mean Square Error: 0.3333333
FA, means and individual ( 95 %) CI
PESO std r LCL UCL Min Max
A1 15.66667 4.226898 6 12.67179 18.66155 10 20
A2 17.66667 5.006662 6 14.67179 20.66155 12 24
Alpha: 0.05 ; DF Error: 1
Critical Value of t: 12.7062
least Significant Difference: 4.235402
Treatments with the same letter are not significantly different.
PESO groups
A2 17.66667 a
A1 15.66667 a
Comparaciones de medias para los efectos principales del Factor B
data %>%with(LSD.test( PESO, # Cambiar según nombre de variable respuesta FB, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dcl), MSerror = dvmisc::get_mse(modelo.dcl),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: PESO ~ FB
LSD t Test for PESO
Mean Square Error: 1.333333
FB, means and individual ( 95 %) CI
PESO std r LCL UCL Min Max
B1 13.75 3.862210 4 12.14702 15.35298 10 18
B2 20.25 3.862210 4 18.64702 21.85298 16 24
B3 16.00 4.082483 4 14.39702 17.60298 12 20
Alpha: 0.05 ; DF Error: 4
Critical Value of t: 2.776445
least Significant Difference: 2.266958
Treatments with the same letter are not significantly different.
PESO groups
B2 20.25 a
B3 16.00 b
B1 13.75 b
Precaución
Si el análisis de varianza arroja que los niveles de un factor son estadísticamente similares entre sí, entonces no es necesario realizar una prueba de comparación de medias y por ende todos estos niveles pertenecen a un mismo grupo de significancia “a”.
Comparaciones de medias para las interacciones
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
B1 vs B3:
\(H_0: \mu_{B1} - \mu_{B3} = 0\)
\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)
B2 vs B3:
\(H_0: \mu_{B2} - \mu_{B3} = 0\)
\(H_1: \mu_{B2} - \mu_{B3} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
Análisis de varianza para interacción de dos factores con el paquete phia
Comparación de los niveles de B dentro de cada nivel de A
Comparaciones de medias de los niveles de B dentro de cada nivel de A
FA_A1 <- data %>%filter(FA=="A1")FA_A2 <- data %>%filter(FA=="A2")fcomp1 <-function(x){ comp <-LSD.test(x$PESO, # Cambiar según nombre de variable respuesta x$FB, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dcl), MSerror = dvmisc::get_mse(modelo.dcl),alpha =0.05,group=TRUE,main =NULL,console=FALSE)return(comp[[5]] %>%rename("Peso"="x$PESO") # Cambiar según nombre de variable respuesta )}
Comparaciones de los niveles de B dentro del nivel A1
fcomp1(FA_A1)
Peso groups
B3 19.5 a
B2 17.0 a
B1 10.5 b
Comparaciones de los niveles de B dentro del nivel A2
fcomp1(FA_A2)
Peso groups
B2 23.5 a
B1 17.0 b
B3 12.5 c
Comparaciones de medias de los niveles de B dentro de cada nivel de A
FB_B1 <- data %>%filter(FB=="B1")FB_B2 <- data %>%filter(FB=="B2")FB_B3 <- data %>%filter(FB=="B3")fcomp2 <-function(x){ comp <-LSD.test(x$PESO, # Cambiar según nombre de variable respuesta x$FA, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dcl), MSerror = dvmisc::get_mse(modelo.dcl),alpha =0.05,group=TRUE,main =NULL,console=FALSE)return(comp[[5]] %>%rename("Peso"="x$PESO") # Cambiar según nombre de variable respuesta )}
Comparaciones de los niveles de A dentro del nivel B1
fcomp2(FB_B1)
Peso groups
A2 17.0 a
A1 10.5 b
Comparaciones de los niveles de A dentro del nivel B2
fcomp2(FB_B2)
Peso groups
A2 23.5 a
A1 17.0 b
Comparaciones de los niveles de A dentro del nivel B3
fcomp2(FB_B3)
Peso groups
A1 19.5 a
A2 12.5 b
Diseño de bloques completos al azar en arreglo de parcelas sub-divididas
concatenar_n_ssplot <-function(vector, n) { s <-NULLfor(i inseq(from=1, to=length(vector), by=n)) { indices <- i:(i+n-1) grupo <- vector[indices] combinacion <-paste(grupo, collapse=" ") s <-c(s, combinacion) }return(s)}q <-concatenar_n_ssplot(vector = book$FB,# n es el número de parcelas medianas por parcela principaln =2)print(t(matrix(q,c(length(FA),r))))
Durbin-Watson test
data: modelo.dbca
DW = 1.7507, p-value = 0.4711
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.98693, p-value = 0.2291
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.50255, p-value = 0.2024
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.073687, p-value = 0.06947
Asymptotic one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.070607, p-value = 0.5114
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.078414, p-value = 0.2154
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 15.111, p-value = 0.2354
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.9844, p-value = 0.1135
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 3.795665, Df = 1, p = 0.051385
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 77.15, df = 54, p-value = 0.02102
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 92.635, df = 54, p-value = 0.000839
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
-----------------------------
DF = 1
Chi2 = 3.795665
Prob > Chi2 = 0.05138547
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
\(H_1: \text{En al menos una interacción entre el factor A y el factor C el } (\tau\gamma) \text{ es diferente a los demás.}\)
\(H_1: (\tau\gamma)_{ik} \neq 0\text{; en al menos una interacción entre el factor A y el factor C.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(A:parcelagrande)” en el caso de efectos aleatorios o “Error(A:bloque)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos.
# A tibble: 8 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 bloque bloque 2 0.732 0.366 NA NA
2 bloque:nitrogeno nitrogeno 4 61.6 15.4 27.7 9.73e- 5
3 bloque:nitrogeno Residuals 8 4.45 0.556 NA NA
4 bloque:nitrogeno:manejo manejo 2 42.9 21.5 94.8 3.40e-13
5 bloque:nitrogeno:manejo Residuals 28 6.34 0.226 NA NA
6 Within variedad 2 206. 103. 221. 2.60e-33
7 Within nitrogeno:v… 8 14.1 1.77 3.79 8.09e- 4
8 Within Residuals 80 37.3 0.466 NA NA
Valor de la tabla de F para el factor A con una significancia de 0.05.
qf(0.95, 4, 8)
[1] 3.837853
Valor de la tabla de F para el factor B con una significancia de 0.05.
qf(0.95, 2, 28)
[1] 3.340386
Valor de la tabla de F para la interacción A:B con una significancia de 0.05.
qf(0.95, 2, 80)
[1] 3.110766
Conclusión.
Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor A tiene un efecto estadísticamente distinto de 0.
Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tiene un efecto estadísticamente distinto de 0.
Con respecto al Factor C: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor C tiene un efecto estadísticamente distinto de 0.
Con respecto a la interacción entre el Factor A y Factor C: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos una interacción entre un nivel del factor A y un nivel del factor C existe un efecto de antagonismo o sinergismo sobre el rendimiento.
Comparaciones de medias para los efectos principales del Factor A
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta nitrogeno, # Cambiar según nombre de variable independienteDFerror =get_df_ea_spspplot(aov.dbca, "nitrogeno"), MSerror =get_mse_ea_spspplot(aov.dbca, "nitrogeno"),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ nitrogeno
LSD t Test for rdto
Mean Square Error: 0.5564188
nitrogeno, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
0 5.384704 1.206963 27 5.053665 5.715743 3.320 8.020
110 6.937370 1.474547 27 6.606331 7.268409 4.246 9.660
140 7.233926 1.968153 27 6.902887 7.564965 3.132 10.360
50 6.220333 1.615861 27 5.889294 6.551372 3.188 9.942
80 6.995741 1.371302 27 6.664702 7.326780 4.422 9.320
Alpha: 0.05 ; DF Error: 8
Critical Value of t: 2.306004
least Significant Difference: 0.4681598
Treatments with the same letter are not significantly different.
rdto groups
140 7.233926 a
80 6.995741 a
110 6.937370 a
50 6.220333 b
0 5.384704 c
Comparaciones de medias para los efectos principales del Factor B
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta manejo, # Cambiar según nombre de variable independienteDFerror =get_df_eb_spspplot(aov.dbca, "nitrogeno", "manejo"), MSerror =get_mse_eb_spspplot(aov.dbca, "nitrogeno", "manejo"),alpha =0.05,group =TRUE,main =NULL,console=TRUE))
Study: rdto ~ manejo
LSD t Test for rdto
Mean Square Error: 0.2264039
manejo, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
m1 5.900378 1.486531 45 5.755082 6.045673 3.132 9.314
m2 6.486156 1.621936 45 6.340860 6.631451 3.625 9.680
m3 7.276711 1.635020 45 7.131416 7.422007 4.225 10.360
Alpha: 0.05 ; DF Error: 28
Critical Value of t: 2.048407
least Significant Difference: 0.2054788
Treatments with the same letter are not significantly different.
rdto groups
m3 7.276711 a
m2 6.486156 b
m1 5.900378 c
Comparaciones de medias para los efectos principales del Factor C
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta variedad, # Cambiar según nombre de variable independienteDFerror =df.residual(modelo.dbca), MSerror = dvmisc::get_mse(modelo.dbca),alpha =0.05,group =TRUE,main =NULL,console=TRUE))
Study: rdto ~ variedad
LSD t Test for rdto
Mean Square Error: 0.4660436
variedad, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 5.126822 0.971530 45 4.924299 5.329345 3.132 7.309
2 6.396133 1.065752 45 6.193611 6.598656 4.166 8.950
3 8.140289 1.314438 45 7.937766 8.342812 5.244 10.360
Alpha: 0.05 ; DF Error: 80
Critical Value of t: 1.990063
least Significant Difference: 0.2864105
Treatments with the same letter are not significantly different.
rdto groups
3 8.140289 a
2 6.396133 b
1 5.126822 c
Precaución
Si el análisis de varianza arroja que los niveles de un factor son estadísticamente similares entre sí, entonces no es necesario realizar una prueba de comparación de medias y por ende todos estos niveles pertenecen a un mismo grupo de significancia “a”.
Comparaciones de medias para las interacciones FAxFB
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
A1 vs A3:
\(H_0: \mu_{A1} - \mu_{A3} = 0\)
\(H_1: \mu_{A1} - \mu_{A3} \neq 0\)
“…”
A4 vs A5:
\(H_0: \mu_{A4} - \mu_{A5} = 0\)
\(H_1: \mu_{A4} - \mu_{A5} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
B1 vs B3:
\(H_0: \mu_{B1} - \mu_{B3} = 0\)
\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)
B2 vs B3:
\(H_0: \mu_{B2} - \mu_{B3} = 0\)
\(H_1: \mu_{B2} - \mu_{B3} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
datasp <- data %>%group_by(bloque, nitrogeno, manejo) %>%summarise(rdto =mean(rdto, na.rm = T)) %>%ungroup()slice_head(datasp, n =5) %>%gt()
Comparaciones de medias de los niveles de B dentro de cada nivel de A
FA_A1 <- data2 %>%filter(Condición=="C0")FA_A2 <- data2 %>%filter(Condición=="C1")fcomp1 <-function(x){ comp <-LSD.test(x$Rend, # Cambiar según nombre de variable respuesta x$Dosis, # Cambiar según nombre de variable independienteDFerror =get_df_eb_spspplot(aov.dbca.ej2, "Condición", "Dosis"), MSerror =get_mse_eb_spspplot(aov.dbca.ej2, "Condición", "Dosis"),alpha =0.05,group=TRUE,main =NULL,console=FALSE)return(comp[[5]] %>%rename("Rend"="x$Rend") # Cambiar según nombre de variable respuesta )}
Comparaciones de medias de los niveles de B dentro de cada nivel de A
fcomp1(FA_A1)
Rend groups
D2 12113.510 a
D3 10733.627 b
D0 9424.306 c
D1 9353.512 c
Comparaciones de los niveles de B dentro del nivel A2
fcomp1(FA_A2)
Rend groups
D2 9228.209 a
D3 8475.458 b
D0 7600.958 c
D1 7326.637 c
$D0
y groups
C0 9424.306 a
C1 7600.958 b
$D1
y groups
C0 9353.512 a
C1 7326.637 b
$D2
y groups
C0 12113.510 a
C1 9228.209 b
$D3
y groups
C0 10733.627 a
C1 8475.458 b
$D0
y groups
C0 9424.306 a
C1 7600.958 b
$D1
y groups
C0 9353.512 a
C1 7326.637 b
$D2
y groups
C0 12113.510 a
C1 9228.209 b
$D3
y groups
C0 10733.627 a
C1 8475.458 b
Comparación de los niveles del factor B dentro de cada nivel del factor A
multcomp.test_2factors <-function(object, respuesta, factor_name1, factor_name2, test, plot, splot, aov){# Función auxiliar para aplicar la prueba de comparaciones múltiples a un data frame multcomp_df <-function(df, respuesta, factor_name, test, plot, splot, aov){if (test =="LSD") { comp <-LSD.test(df[[respuesta]], df[[factor_name]],DFerror =get_df_eb_spspplot(aov, plot, splot), MSerror =get_mse_eb_spspplot(aov, plot, splot),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="HSD") { comp <-HSD.test(df[[respuesta]], df[[factor_name]],DFerror =get_df_eb_spspplot(aov, plot, splot), MSerror =get_mse_eb_spspplot(aov, plot, splot),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="duncan") { comp <-duncan.test(df[[respuesta]], df[[factor_name]],DFerror =get_df_eb_spspplot(aov, plot, splot), MSerror =get_mse_eb_spspplot(aov, plot, splot),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } elseif (test =="SNK") { comp <-SNK.test(df[[respuesta]], df[[factor_name]],DFerror =get_df_eb_spspplot(aov, plot, splot), MSerror =get_mse_eb_spspplot(aov, plot, splot),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } else {stop("Test no válido. Los tests disponibles son LSD, HSD, duncan, y SNK.") }return(comp %>%rename("y"="df[[respuesta]]") ) }# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters1 comp_filters1 <-lapply(object[[1]], function(df){multcomp_df(df, respuesta, factor_name2, test, plot, splot, aov) %>%arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) })# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters2 comp_filters2 <-lapply(object[[2]], function(df){multcomp_df(df, respuesta, factor_name1, test, plot, splot, aov) %>% dplyr::arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) %>% dplyr::mutate(groups =toupper(groups)) })row.names(comp_filters1) <-NULLrow.names(comp_filters2) <-NULL# Retornar una lista con las comparaciones múltiples para cada data frame result <-list() result[[as.name(substitute(factor_name1))]] <- comp_filters1 result[[as.name(substitute(factor_name2))]] <- comp_filters2return(#list(#"Comparación de los niveles del factor B dentro de cada nivel del factor A", result#[[1]],# "Comparación de los niveles del factor A dentro de cada nivel del factor B",# result[[2]]# ) )}
$Condición
$Condición$C0
x y groups
1 D0 9424.306 c
2 D1 9353.512 c
3 D2 12113.510 a
4 D3 10733.627 b
$Condición$C1
x y groups
1 D0 7600.958 c
2 D1 7326.637 c
3 D2 9228.209 a
4 D3 8475.458 b
$Dosis
$Dosis$D0
x y groups
1 C0 9424.306 A
2 C1 7600.958 B
$Dosis$D1
x y groups
1 C0 9353.512 A
2 C1 7326.637 B
$Dosis$D2
x y groups
1 C0 12113.510 A
2 C1 9228.209 B
$Dosis$D3
x y groups
1 C0 10733.627 A
2 C1 8475.458 B
# Convertir la lista a un data framedf <- result.comp %>% reshape2::melt() %>%ungroup() %>% dplyr::select(-variable) %>%relocate("Factor"="L1","Nivel"="L2", x,"y"="value", groups)df
Factor Nivel x y groups
1 Condición C0 D0 9424.306 c
2 Condición C0 D1 9353.512 c
3 Condición C0 D2 12113.510 a
4 Condición C0 D3 10733.627 b
5 Condición C1 D0 7600.958 c
6 Condición C1 D1 7326.637 c
7 Condición C1 D2 9228.209 a
8 Condición C1 D3 8475.458 b
9 Dosis D0 C0 9424.306 A
10 Dosis D0 C1 7600.958 B
11 Dosis D1 C0 9353.512 A
12 Dosis D1 C1 7326.637 B
13 Dosis D2 C0 12113.510 A
14 Dosis D2 C1 9228.209 B
15 Dosis D3 C0 10733.627 A
16 Dosis D3 C1 8475.458 B
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(groups =paste0(groups.y,groups.x)) %>% dplyr::select(!!level1,!!level2, y, groups)# devolver una lista con los dos subconjuntosreturn(df)}
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(y =paste0(round(y,2), " ", groups.y, groups.x)) %>% dplyr::select(!!level1,!!level2, y)# devolver una lista con los dos subconjuntosreturn(df)}
distribuir_vector_matriz <-function(vector, niveles_FA, niveles_FB, num_bloques) {# Verificar que el número de niveles de FB sea menor o igual a la mitad del número de bloquesstopifnot(niveles_FB <= num_bloques/2)# Calcular el número de columnas por bloque columnas_bloque <- niveles_FA# Calcular el número de filas por bloque filas_bloque <- niveles_FB# Calcular el número total de filas num_filas <- num_bloques * filas_bloque# Calcular el número de bloques completos num_bloques_completos <- num_filas %/% filas_bloque# Crear una matriz vacía matriz <-matrix(NA, nrow = num_filas, ncol = columnas_bloque)# Distribuir los valores en la matrizfor (i in1:num_bloques_completos) { inicio_fila <- (i -1) * filas_bloque +1 fin_fila <- i * filas_bloque inicio_vector <- (i -1) * columnas_bloque * filas_bloque +1 fin_vector <- i * columnas_bloque * filas_bloque matriz[inicio_fila:fin_fila, ] <-matrix(vector[inicio_vector:fin_vector], nrow = filas_bloque, byrow =TRUE) %>%t() }# Manejar los bloques incompletosif (num_filas %% filas_bloque !=0) { inicio_fila <- num_bloques_completos * filas_bloque +1 fin_fila <- num_filas inicio_vector <- num_bloques_completos * columnas_bloque * filas_bloque +1 fin_vector <-length(vector) matriz[inicio_fila:fin_fila, ] <-matrix(vector[inicio_vector:fin_vector], nrow = fin_fila - inicio_fila +1, byrow =TRUE) }return(matriz)}
distribuir_vector_matriz <-function(vector, niveles_FA, niveles_FB, num_bloques) {# Verificar que el número de niveles de FB sea menor o igual a la mitad del número de bloquesstopifnot(niveles_FB <= num_bloques/2)# Calcular el número de columnas por bloque columnas_bloque <- niveles_FA# Calcular el número de filas por bloque filas_bloque <- niveles_FB# Calcular el número total de filas num_filas <- filas_bloque# Calcular el número total de filas num_columnas <- columnas_bloque * num_bloques# Calcular el número de bloques completos num_bloques_completos <- num_filas %/% filas_bloque# Crear una matriz vacía matriz <-matrix(NA, nrow = num_filas, ncol = num_columnas)# Distribuir los valores en la matrizfor (i in1:num_bloques_completos) { inicio_fila <- (i -1) * filas_bloque +1 fin_fila <- i * filas_bloque inicio_vector <- (i -1) * columnas_bloque * filas_bloque +1 fin_vector <- i * columnas_bloque * filas_bloque matriz[inicio_fila:fin_fila, ] <-matrix(vector[inicio_vector:fin_vector], nrow = filas_bloque, byrow =TRUE) %>%t() }# Manejar los bloques incompletosif (num_filas %% filas_bloque !=0) { inicio_fila <- num_bloques_completos * filas_bloque +1 fin_fila <- num_filas inicio_vector <- num_bloques_completos * columnas_bloque * filas_bloque +1 fin_vector <-length(vector) matriz[inicio_fila:fin_fila, ] <-matrix(vector[inicio_vector:fin_vector], nrow = fin_fila - inicio_fila +1, byrow =TRUE) }return(matriz)}
Durbin-Watson test
data: modelo.dbca
DW = 2.6957, p-value = 0.278
alternative hypothesis: true autocorrelation is not 0
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los residuos del rendimiento son completamente aleatorios e independientes.
Normalidad de residuos
\(H_0: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
\(H_1: \text{La distribución de los residuos del rendimiento es similar a la función normal}\)
shapiro.test(rstudent(modelo.dbca))
Shapiro-Wilk normality test
data: rstudent(modelo.dbca)
W = 0.92973, p-value = 0.09615
ad.test(rstudent(modelo.dbca))
Anderson-Darling normality test
data: rstudent(modelo.dbca)
A = 0.73842, p-value = 0.04692
lillie.test(rstudent(modelo.dbca))
Lilliefors (Kolmogorov-Smirnov) normality test
data: rstudent(modelo.dbca)
D = 0.14393, p-value = 0.2231
Exact one-sample Kolmogorov-Smirnov test
data: rstudent(modelo.dbca)
D = 0.096589, p-value = 0.9627
alternative hypothesis: two-sided
cvm.test(rstudent(modelo.dbca))
Cramer-von Mises normality test
data: rstudent(modelo.dbca)
W = 0.12618, p-value = 0.04565
pearson.test(rstudent(modelo.dbca))
Pearson chi-square normality test
data: rstudent(modelo.dbca)
P = 8, p-value = 0.1562
sf.test(rstudent(modelo.dbca))
Shapiro-Francia normality test
data: rstudent(modelo.dbca)
W = 0.90414, p-value = 0.02824
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la distribución de los residuos del rendimiento es similar a la función normal o gaussiana.
Homocedasticidad
\(H_0\): La varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento
\(H_1\): La varianza del rendimiento no es constante con respecto a los valores ajustados del rendimiento
ncvTest(modelo.dbca)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 2.73003, Df = 1, p = 0.098477
bptest(modelo.dbca)
studentized Breusch-Pagan test
data: modelo.dbca
BP = 24, df = 17, p-value = 0.1194
bptest(modelo.dbca, studentize = F)
Breusch-Pagan test
data: modelo.dbca
BP = 21.652, df = 17, p-value = 0.1985
olsrr::ols_test_breusch_pagan(modelo.dbca)
Breusch Pagan Test for Heteroskedasticity
-----------------------------------------
Ho: the variance is constant
Ha: the variance is not constant
Data
--------------------------------
Response : rdto
Variables: fitted values of rdto
Test Summary
-----------------------------
DF = 1
Chi2 = 2.73003
Prob > Chi2 = 0.09847741
Conclusión. A un nivel de significancia de 0.1, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, la varianza del rendimiento es constante con respecto a los valores ajustados del rendimiento.
Recomendación. Debido a que se cumple con el supuesto de homocedasticidad, para evaluar los efectos de los tratamientos con respecto al rendimiento, se debe proceder a realizar el análisis de varianza.
Estadísticas globales
modelo.dbca %>%gvlma()
Call:
lm(formula = rdto ~ trat * aplic + bloque/trat + bloque/aplic +
bloque, data = data)
Coefficients:
(Intercept) trat2 trat3 aplic2 aplic3
9.208e+00 -2.250e+00 -1.375e+00 4.500e+00 -3.333e+00
aplic4 bloque2 trat2:aplic2 trat3:aplic2 trat2:aplic3
-5.000e+00 -1.417e+00 2.000e+00 -5.000e-01 1.000e+00
trat3:aplic3 trat2:aplic4 trat3:aplic4 trat2:bloque2 trat3:bloque2
-2.584e-15 3.000e+00 3.000e+00 1.500e+00 7.500e-01
aplic2:bloque2 aplic3:bloque2 aplic4:bloque2
3.000e+00 3.667e+00 -1.000e+00
ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
Level of Significance = 0.05
Call:
gvlma(x = .)
Value p-value Decision
Global Stat 5.146e+00 0.27259 Assumptions acceptable.
Skewness 1.021e-31 1.00000 Assumptions acceptable.
Kurtosis 3.828e-02 0.84488 Assumptions acceptable.
Link Function 4.376e+00 0.03645 Assumptions NOT satisfied!
Heteroscedasticity 7.324e-01 0.39210 Assumptions acceptable.
\(H_1: \text{En al menos una interacción entre el factor A y el factor B el } (\tau\beta) \text{ es diferente a los demás.}\)
\(H_1: (\tau\beta)_{ij} \neq 0\text{; en al menos una interacción entre el factor A y el factor B.}\)
Precaución
Si se ignora el diseño experimental, se obtienen los siguientes resultados, incorrectos. Observe que los grados de libertad del error es mayor, lo que facilita la detección de diferencias que en realidad no existen.
Necesitamos especificar correctamente el término de error para el factor A. Debe tenerse en cuenta que la forma del Error es “Error(bloque:A)” en el caso de efectos aleatorios o “Error(bloque/A)” para el caso de efectos fijos y puede cambiar dependiendo de la disposición de los datos. La clave es conocer los grados de libertad correcto para saber que se obtienen los resultados correctos. En el factor B se trabaja con “Error(bloque:B)” en el caso de efectos aleatorios o “Error(bloque/B)” para el caso de efectos fijos. Para la interacción se trabaja con “Error(bloque/A/B)” para el caso de efectos fijos o con “Error(bloque:A:B)” en el caso de efectos aleatorios.
aov(rdto ~ bloque + trat * aplic +Error(bloque/trat+bloque/aplic), data = data) -> aov.dbcasummary(aov.dbca)
Error: bloque
Df Sum Sq Mean Sq
bloque 1 3.375 3.375
Error: bloque:trat
Df Sum Sq Mean Sq F value Pr(>F)
trat 2 0.75 0.375 0.333 0.75
Residuals 2 2.25 1.125
Error: bloque:aplic
Df Sum Sq Mean Sq F value Pr(>F)
aplic 3 330.1 110.04 14.28 0.0279 *
Residuals 3 23.1 7.71
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
trat:aplic 6 11.25 1.8750 1.957 0.217
Residuals 6 5.75 0.9583
data %>%with(agricolae::strip.plot(BLOCK = bloque,COL = trat,ROW = aplic,Y = rdto))
ANALYSIS STRIP PLOT: rdto
Class level information
trat : 1 2 3
aplic : 3 1 4 2
bloque : 1 2
Number of observations: 24
model Y: rdto ~ bloque + trat + Ea + aplic + Eb + aplic:trat + Ec
Analysis of Variance Table
Response: rdto
Df Sum Sq Mean Sq F value Pr(>F)
bloque 1 3.37 3.375
trat 2 0.75 0.375 0.3333 0.75000
Ea 2 2.25 1.125
aplic 3 330.13 110.042 14.2757 0.02787 *
Eb 3 23.12 7.708
aplic:trat 6 11.25 1.875 1.9565 0.21719
Ec 6 5.75 0.958
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
cv(a) = 13.1 %, cv(b) = 34.2 %, cv(c) = 12 %, Mean = 8.125
broom::tidy(aov.dbca)
# A tibble: 7 × 7
stratum term df sumsq meansq statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 bloque bloque 1 3.37 3.37 NA NA
2 bloque:trat trat 2 0.750 0.375 0.333 0.75
3 bloque:trat Residuals 2 2.25 1.13 NA NA
4 bloque:aplic aplic 3 330. 110. 14.3 0.0279
5 bloque:aplic Residuals 3 23.1 7.71 NA NA
6 Within trat:aplic 6 11.3 1.88 1.96 0.217
7 Within Residuals 6 5.75 0.958 NA NA
Valor de la tabla de F para el factor A con una significancia de 0.05.
qf(0.95, 2, 2)
[1] 19
Valor de la tabla de F para el factor B con una significancia de 0.05.
qf(0.95, 3, 3)
[1] 9.276628
Valor de la tabla de F para la interacción A:B con una significancia de 0.05.
qf(0.95, 6, 6)
[1] 4.283866
Conclusión.
Con respecto al Factor A: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, los niveles del factor A tiene un efecto estadísticamente similar de 0.
Con respecto al Factor B: A un nivel de significancia de 0.05, se concluye que existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, al menos un nivel del factor B tiene un efecto estadísticamente distinto de 0.
Con respecto a la interacción entre el Factor A y Factor B: A un nivel de significancia de 0.05, se concluye que no existe suficiente evidencia estadística para rechazar la hipótesis nula, por lo tanto, no existió un efecto de antagonismo o sinergismo en alguna interacción dada por los niveles de A y de B.
Comparaciones de medias para los efectos principales del Factor A
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta trat, # Cambiar según nombre de variable independienteDFerror =get_dfe_stplot(aov.dbca, "trat"), MSerror =get_mse_stplot(aov.dbca, "trat"),alpha =0.05,group=TRUE,main =NULL,console=TRUE))
Study: rdto ~ trat
LSD t Test for rdto
Mean Square Error: 1.125
trat, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 8.250 4.496030 8 6.636505 9.863495 2 15
2 8.250 4.682795 8 6.636505 9.863495 3 17
3 7.875 3.399054 8 6.261505 9.488495 5 14
Alpha: 0.05 ; DF Error: 2
Critical Value of t: 4.302653
least Significant Difference: 2.281826
Treatments with the same letter are not significantly different.
rdto groups
1 8.250 a
2 8.250 a
3 7.875 a
Comparaciones de medias para los efectos principales del Factor B
data %>%with(LSD.test( rdto, # Cambiar según nombre de variable respuesta aplic, # Cambiar según nombre de variable independienteDFerror =get_dfe_stplot(aov.dbca, "aplic"), MSerror =get_mse_stplot(aov.dbca, "aplic"),alpha =0.05,group =TRUE,main =NULL,console=TRUE))
Study: rdto ~ aplic
LSD t Test for rdto
Mean Square Error: 7.708333
aplic, means and individual ( 95 %) CI
rdto std r LCL UCL Min Max
1 7.666667 0.8164966 6 4.0595042 11.273829 7 9
2 14.166667 1.7224014 6 10.5595042 17.773829 12 17
3 6.500000 1.8708287 6 2.8928375 10.107162 4 9
4 4.166667 1.4719601 6 0.5595042 7.773829 2 6
Alpha: 0.05 ; DF Error: 3
Critical Value of t: 3.182446
least Significant Difference: 5.101298
Treatments with the same letter are not significantly different.
rdto groups
2 14.166667 a
1 7.666667 b
3 6.500000 b
4 4.166667 b
Precaución
Si el análisis de varianza arroja que los niveles de un factor son estadísticamente similares entre sí, entonces no es necesario realizar una prueba de comparación de medias y por ende todos estos niveles pertenecen a un mismo grupo de significancia “a”.
Comparaciones de medias para las interacciones FAxFB
Para los niveles del factor A dentro del nivel B1:
A1 vs A2:
\(H_0: \mu_{A1} - \mu_{A2} = 0\)
\(H_1: \mu_{A1} - \mu_{A2} \neq 0\)
A1 vs A3:
\(H_0: \mu_{A1} - \mu_{A3} = 0\)
\(H_1: \mu_{A1} - \mu_{A3} \neq 0\)
A2 vs A3:
\(H_0: \mu_{A4} - \mu_{A5} = 0\)
\(H_1: \mu_{A4} - \mu_{A5} \neq 0\)
Para los niveles del factor B dentro del nivel A1:
B1 vs B2:
\(H_0: \mu_{B1} - \mu_{B2} = 0\)
\(H_1: \mu_{B1} - \mu_{B2} \neq 0\)
B1 vs B3:
\(H_0: \mu_{B1} - \mu_{B3} = 0\)
\(H_1: \mu_{B1} - \mu_{B3} \neq 0\)
B1 vs B4:
\(H_0: \mu_{B1} - \mu_{B4} = 0\)
\(H_1: \mu_{B1} - \mu_{B4} \neq 0\)
(…)
B3 vs B4:
\(H_0: \mu_{B3} - \mu_{B4} = 0\)
\(H_1: \mu_{B3} - \mu_{B4} \neq 0\)
NOTA: Repetir este proceso para cada nivel de A y cada nivel de B.
multcomp.test_2factors <-function(object, respuesta, factor_name1, factor_name2, test, aov){# Función auxiliar para aplicar la prueba de comparaciones múltiples a un data frame multcomp_df <-function(df, respuesta, factor_name, test, aov){if (test =="LSD") { comp <-LSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="HSD") { comp <-HSD.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[5]] } elseif (test =="duncan") { comp <-duncan.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } elseif (test =="SNK") { comp <-SNK.test(df[[respuesta]], df[[factor_name]],DFerror =df.residual(aov), MSerror = dvmisc::get_mse(aov),alpha =0.05,group=TRUE,main =NULL,console=FALSE)[[6]] } else {stop("Test no válido. Los tests disponibles son LSD, HSD, duncan, y SNK.") }return(comp %>%rename("y"="df[[respuesta]]") ) }# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters1 comp_filters1 <-lapply(object[[1]], function(df){multcomp_df(df, respuesta, factor_name2, test, aov) %>%arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) })# Aplicar la prueba de comparaciones múltiples a todos los data frames dentro de filters2 comp_filters2 <-lapply(object[[2]], function(df){multcomp_df(df, respuesta, factor_name1, test, aov) %>% dplyr::arrange(row.names(.)) %>%rownames_to_column(var ="x") %>%relocate(x) %>% dplyr::mutate(groups =toupper(groups)) })row.names(comp_filters1) <-NULLrow.names(comp_filters2) <-NULL# Retornar una lista con las comparaciones múltiples para cada data frame result <-list() result[[as.name(substitute(factor_name1))]] <- comp_filters1 result[[as.name(substitute(factor_name2))]] <- comp_filters2return(#list(#"Comparación de los niveles del factor B dentro de cada nivel del factor A", result#[[1]],# "Comparación de los niveles del factor A dentro de cada nivel del factor B",# result[[2]]# ) )}
$trat
$trat$`1`
x y groups
1 1 8.5 b
2 2 14.5 a
3 3 7.0 b
4 4 3.0 c
$trat$`2`
x y groups
1 1 7.0 b
2 2 15.0 a
3 3 6.5 bc
4 4 4.5 c
$trat$`3`
x y groups
1 1 7.5 b
2 2 13.0 a
3 3 6.0 bc
4 4 5.0 c
$aplic
$aplic$`1`
x y groups
1 1 8.5 A
2 2 7.0 A
3 3 7.5 A
$aplic$`2`
x y groups
1 1 14.5 A
2 2 15.0 A
3 3 13.0 A
$aplic$`3`
x y groups
1 1 7.0 A
2 2 6.5 A
3 3 6.0 A
$aplic$`4`
x y groups
1 1 3.0 A
2 2 4.5 A
3 3 5.0 A
# Convertir la lista a un data framedf <- result.comp %>% reshape2::melt() %>%ungroup() %>% dplyr::select(-variable) %>%relocate("Factor"="L1","Nivel"="L2", x,"y"="value", groups)df
Factor Nivel x y groups
1 trat 1 1 8.5 b
2 trat 1 2 14.5 a
3 trat 1 3 7.0 b
4 trat 1 4 3.0 c
5 trat 2 1 7.0 b
6 trat 2 2 15.0 a
7 trat 2 3 6.5 bc
8 trat 2 4 4.5 c
9 trat 3 1 7.5 b
10 trat 3 2 13.0 a
11 trat 3 3 6.0 bc
12 trat 3 4 5.0 c
13 aplic 1 1 8.5 A
14 aplic 1 2 7.0 A
15 aplic 1 3 7.5 A
16 aplic 2 1 14.5 A
17 aplic 2 2 15.0 A
18 aplic 2 3 13.0 A
19 aplic 3 1 7.0 A
20 aplic 3 2 6.5 A
21 aplic 3 3 6.0 A
22 aplic 4 1 3.0 A
23 aplic 4 2 4.5 A
24 aplic 4 3 5.0 A
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# # Convertir level1 y level2 en nombres simbólicos# level1 <- as.name(level1)# level2 <- as.name(level2)# # filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1, level2)) %>% dplyr::mutate(groups =paste0(groups.y,groups.x)) %>% dplyr::select(!!level1,!!level2, y, groups) #%>%# rename(!!level1 := x,# !!level2 := Nivel)# devolver una lista con los dos subconjuntosreturn(df)}
# función que filtra por los dos niveles especificados por el usuario y devuelve dos subconjuntos de datoscreate_report <-function(df, level1, level2) {# filtrar por el primer nivel subset1 <- df %>%filter(Factor == level1) %>%rename(!!level2 := x,!!level1 := Nivel) %>% dplyr::select(-c(Factor))# filtrar por el segundo nivel subset2 <- df %>%filter(Factor == level2) %>%rename(!!level1 := x,!!level2 := Nivel) %>% dplyr::select(-c(Factor,y)) df <- subset1 %>% dplyr::left_join(subset2,by =c(level1,level2)) %>% dplyr::mutate(y =paste0(round(y,2), " ", groups.y, groups.x)) %>% dplyr::select(!!level1,!!level2, y)# devolver una lista con los dos subconjuntosreturn(df)}
------------------------------------------------------------------------
Legend:
FACTOR 1 (Whole plot): Tratamiento
FACTOR 2 (strip-plot): Aplicación
------------------------------------------------------------------------
------------------------------------------------------------------------
Analysis of variance table
------------------------------------------------------------------------
DF SS MS Fc Pr(>Fc)
Block 1 3.38 3.375 0.4286 0.57992
Tratamiento 2 0.75 0.375 0.3333 0.75000
Error a 2 2.25 1.125
Aplicación 3 330.12 110.042 14.2757 0.02787 *
Error b 3 23.12 7.708
Tratamiento*Aplicación 6 11.25 1.875 1.9565 0.21719
Error c 6 5.75 0.958
Total 23 376.62
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
------------------------------------------------------------------------
CV 1 = 13.05428 %
CV 2 = 34.17094 %
CV 3 = 12.04855 %
------------------------------------------------------------------------
Shapiro-Wilk normality test (Error b)
p-value: 0.9973946
According to Shapiro-Wilk normality test at 5% of significance, residuals can be considered normal.
------------------------------------------------------------------------
No significant interaction: analyzing the main effects
------------------------------------------------------------------------
Tratamiento
According to F test, the means of this factor are not different.
------------------------------------------------------------------------
Levels Means
1 1 8.250
2 2 8.250
3 3 7.875
------------------------------------------------------------------------
Aplicación
Duncan's test
------------------------------------------------------------------------
Groups Treatments Means
a 2 14.16667
b 1 7.666667
b 3 6.5
b 4 4.166667
------------------------------------------------------------------------
